Visualize geographic map¶

Created by

¶ Li, Chaonan (李超男) / licn@mtc.edu.cn / Ecological Security and Protection Key Laboratory of Sichuan Province, Mianyang Normal University (绵阳师范学院生态安全与保护四川省重点实验室)

¶ Liu Chi (刘驰) / liuchi0426@126.com / College of Resources and Environment, Fujian Agriculture and Forestry University (福建农林大学资源与环境学院)

¶ Liao, Haijun (廖海君) / lihj@mtc.edu.cn / Engineering Research Center of Chuanxibei RHS Construction at Mianyang Normal University of Sichuan Province (绵阳师范学院川西北乡村人居环境建设工程研究中心)

Reviewed by

¶ Li, Xiangzhen (李香真) / lixz@fafu.edu.cn / College of Resources and Environment, Fujian Agriculture and Forestry University (福建农林大学资源与环境学院)

One of our purposes in developing the microgeo R package is to rapidly visualize the biogeographic patterns of microbial traits onto maps. Hence, we have designed three visualization functions, and all of them return a ggplot2 object. The function of plot_bmap() is able to accept a microgeo-compatible SpatialPolygonsDataFrame returned by read_aliyun_map(), trans_map_fmt(), grid_map(), merge_dfs_to_map() or merge_mtx_to_map(), and it only visualize the data included in the SpatialPolygonsDataFrame itself. The plot_nmap() and plot_imap() are designed to visualize the results of nearest neighbour and inverse distance weighting interpolations, respectively.

To add more microbial traits onto maps, we have implemented four functions including add_label(), add_point(), add_spatraster() and add_contour. The add_label() can add any characters or numbers in the native polygons of map or the grids created by grid_map(). For example, we can apply the gradient colors to represent the average shannon index in each grid, and simultaneously add the standard errors into each grid by using the add_label(). The add_point() can add points based on a numeric vector onto a map. By applying such a function, we also can use the size of point to represent the standard errors of shannon index. Meanwhile, we also can use the gradient colors and/or size of point as a proxy of any microbial trait at each sampling site on maps.

In the microgeo R package, the functions responsible for predictions based on machine learning model and some spatial interpolations return a SpatRaster object. To visualize such a type of data onto maps, we have designed a function of add_spatraster(). To convert the data in the SpatRaster to contours that can be further rendered on maps, we have designed a function of add_contours().

Now, let's go through each of these functions and see how to visualize microbial/environmental traits onto maps.

Load required R packages¶

In [1]:
# Load three required packages 
suppressMessages(require("magrittr")) 
require("ggplot2")  %>% suppressMessages()
require("microgeo") %>% suppressMessages()
prev_locale <- Sys.setlocale("LC_CTYPE", "C.UTF-8") # Set locale to `UTF-8` to show Chinese characters in tables

Plot basic map¶

The microgeo R package provides two test datasets for the plot_bmap() function, where the common.map.mean is a SpatialPolygonsDataFrame based on administrative divisions and the gridded.map.mean is a SpatialPolygonsDataFrame processed by grid_map(). Below is a simple example of the plot_bmap() function.

Firstly, Let's check the content of common.map.mean, a microgeo-compatible SpatialPolygonsDataFrame.

In [2]:
# This is a microgeo-compatible SpatialPolygonsDataFrame
data(common.map.mean)
head(common.map.mean@data[,1:12])
A data.frame: 3 × 12
TYPEFMTSNAMEX.CENTERY.CENTERobserved_meanshannon_meanobserved_sdshannon_sdobserved_seshannon_sesample.num
<chr><chr><chr><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><int>
1DataV.GeoAtlasmicrogeo西藏自治区 88.3882831.56375663.84165.846949239.04070.478179210.442540.02088936524
2DataV.GeoAtlasmicrogeo青海省 96.0435335.72640647.70365.837889246.25420.539671211.713120.02566954442
3DataV.GeoAtlasmicrogeo四川省 102.6934530.67454706.84175.974116198.00670.365985911.875650.02195037278
In [3]:
# Change the name of Polygons
common.map.mean$NAME <- c("Tibet", "Qinghai", "Sichuan") 
head(common.map.mean@data[,1:12]) %>% suppressMessages()

Loading required package: sp
A data.frame: 3 × 12
TYPEFMTSNAMEX.CENTERY.CENTERobserved_meanshannon_meanobserved_sdshannon_sdobserved_seshannon_sesample.num
<chr><chr><chr><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><int>
1DataV.GeoAtlasmicrogeoTibet 88.3882831.56375663.84165.846949239.04070.478179210.442540.02088936524
2DataV.GeoAtlasmicrogeoQinghai 96.0435335.72640647.70365.837889246.25420.539671211.713120.02566954442
3DataV.GeoAtlasmicrogeoSichuan102.6934530.67454706.84175.974116198.00670.365985911.875650.02195037278

Then, we use plot_bmap to visualize the common.map.mean.

In [4]:
# Plot a basic map without any fill color 
options(repr.plot.width = 15.4, repr.plot.height = 8.02)
(p1 <- common.map.mean %>% plot_bmap(legend.position = c(0.06, 0.25)) %>% suppressMessages())
No description has been provided for this image
In [5]:
# Plot a basic map with the naive polygons filled by the automatically generated colors based on a character variable
options(repr.plot.width = 15.4, repr.plot.height = 8.02)
(p2 <- common.map.mean %>% plot_bmap(var = 'NAME', fill = 'auto', legend.position = c(0.06, 0.25)))
No description has been provided for this image
In [6]:
# Plot a basic map with the naive polygons filled by the automatically generated colors based on a numeric variable
options(repr.plot.width = 15.4, repr.plot.height = 8.02)
(p3 <- common.map.mean %>% plot_bmap(var = 'shannon_mean', fill = 'auto', legend.position = c(0.06, 0.255)))
No description has been provided for this image
In [7]:
# Plot a basic map with the naive polygons filled by manually set colors based on a character variable
# The display order of polygon names are also manually set
options(repr.plot.width = 15.4, repr.plot.height = 8.02)
(p4 <- common.map.mean %>% plot_bmap(var = 'NAME', fill = RColorBrewer::brewer.pal(8, "Set2")[1:3], 
                                     ord = c('Tibet', 'Qinghai', 'Sichuan'), legend.position = c(0.06, 0.25)))
No description has been provided for this image
In [8]:
# Plot a basic map with the naive polygons filled by manually set colors based on a numeric variable
options(repr.plot.width = 15.4, repr.plot.height = 8.02)
(p5 <- common.map.mean %>% plot_bmap(var = 'shannon_mean', fill = colorRampPalette(RColorBrewer::brewer.pal(11, "BrBG"))(100), 
                                     legend.position = c(0.06, 0.255)))
No description has been provided for this image

Now, let's check the content of gridded.map.mean, a microgeo-compatible SpatialPolygonsDataFrame. The NA means that there are no samples in the grid. The NAME means the grid ID.

In [9]:
# This is a microgeo-compatible SpatialPolygonsDataFrame
data(gridded.map.mean) 
head(gridded.map.mean@data[,1:12])
A data.frame: 6 × 12
TYPEFMTSNAMEX.CENTERY.CENTERobserved_meanshannon_meanobserved_sdshannon_sdobserved_seshannon_sesample.num
<chr><chr><chr><dbl><dbl><dbl><dbl><dbl><dbl><dbl><dbl><int>
1Gridded.Mapmicrogeo1 92.7270827.05839NANANANANANANA
2Gridded.Mapmicrogeo2 93.6024027.08612NANANANANANANA
3Gridded.Mapmicrogeo3102.6813126.73440NANANANANANANA
4Gridded.Mapmicrogeo4101.8870626.72948NANANANANANANA
5Gridded.Mapmicrogeo5 94.0810628.05720NANANANANANANA
6Gridded.Mapmicrogeo6100.4811228.28070NANANANANANANA

Then, we use plot_bmap to visualize the gridded.map.mean.

In [10]:
# Plot a gridded map without any fill color 
options(repr.plot.width = 15.4, repr.plot.height = 8.02)
(p6 <- gridded.map.mean %>% plot_bmap(legend.position = c(0.06, 0.25))) %>% suppressMessages()
No description has been provided for this image
In [11]:
# Plot a gridded map with the grids filled by the automatically generated colors based on a numeric variable
options(repr.plot.width = 15.4, repr.plot.height = 8.02)
(p7 <- gridded.map.mean %>% plot_bmap(var = 'shannon_mean', fill = 'auto', legend.position = c(0.06, 0.255)))
No description has been provided for this image
In [12]:
# Plot a gridded map with the grids filled by manually set colors based on a numeric variable
options(repr.plot.width = 15.4, repr.plot.height = 8.02)
(p8 <- gridded.map.mean %>% plot_bmap(var = 'shannon_mean', fill = rev(colorRampPalette(RColorBrewer::brewer.pal(11, "BrBG"))(100)),
                                      legend.position = c(0.06, 0.255)))
No description has been provided for this image

Plot a nearest neighbour interpolation map¶

We have designed a function of plot_nmap() to visualize the results of nearest neighbour interpolation. Here, we only give a simple example for the visualization. See details in the tutorial of nearest neighbour interpolation.

Firstly, let's check the content of nen.interp.

In [13]:
# Load test data
data(nen.interp)
In [14]:
# The `nen.interp` is a list
class(nen.interp)
'list'
In [15]:
# There are two datasets in the `nen.interp`
names(nen.interp)
  1. 'shannon'
  2. 'observed'

Then, we use plot_nmap to visualize the nen.interp.

In [16]:
# Plot a map for the results of nearest neighbour interpolation [alpha diversity ==> observed] 
options(repr.plot.width = 15.4, repr.plot.height = 8.02)
(p9 <- nen.interp$observed %>% plot_nmap(legend.position = c(0.06, 0.25)))
No description has been provided for this image
In [17]:
# Plot a map for the results of nearest neighbour interpolation [alpha diversity ==> shannon]; Use manually set colors
options(repr.plot.width = 15.4, repr.plot.height = 8.02)
(p10 <- nen.interp$shannon %>% plot_nmap(fill = colorRampPalette(RColorBrewer::brewer.pal(11, "BrBG"))(100), legend.position = c(0.06, 0.255)))
No description has been provided for this image

Plot an inverse distance weighting (IDW) interpolation map¶

The microgeo provides two types of inverse distance weighting (IDW) interpolation methods, one of which is hexagonal interpolation. It generates a sf data.frame and requires a separate plotting function for visualization. To address this issue, we have designed a function of plot_imap(). Here, we only give a simple example for visualization. See details in the tutorial of IDW interpolation.

Firstly, let's check the content of idw.interp.hex.

In [18]:
# Load test data
data(idw.interp.hex)
In [19]:
# The `idw.interp.hex` is a list
class(idw.interp.hex)
'list'
In [20]:
# There are two datasets in the `idw.interp.hex`
names(idw.interp.hex)
  1. 'observed'
  2. 'shannon'

Then, we use plot_imap to visualize the idw.interp.hex.

In [21]:
# Plot a map for the results of inverse distance weighting (IDW) interpolation [alpha diversity ==> observed] 
options(repr.plot.width = 15.4, repr.plot.height = 8.02)
(p11 <- idw.interp.hex$observed %>% plot_imap(legend.position = c(0.06, 0.255)))
No description has been provided for this image
In [22]:
# Plot a map for the results of inverse distance weighting (IDW) interpolation [alpha diversity ==> observed]; Use manually set colors
options(repr.plot.width = 15.4, repr.plot.height = 8.02)
(p12 <- idw.interp.hex$observed %>% plot_imap(fill = colorRampPalette(RColorBrewer::brewer.pal(11, "BrBG"))(100), legend.position = c(0.06, 0.255)))
No description has been provided for this image

Add labels onto a map¶

If you want to add some labels onto maps, the add_label() function may be very useful. With the add_label(), you can add the names of administrative region, grid names, and any character or numeric variables included in SpatialPolygonsDataFrame onto maps. You can even use this function multiple times to add multiple labels onto maps. Here are a few simple examples.

In [23]:
# Add grid NAME on the map
options(repr.plot.width = 15.4 * 2, repr.plot.height = 8.02)
p13 <- p6 %>% add_label(dat = gridded.map.mean@data, lab.var = 'NAME', lon.var = 'X.CENTER', 
                        lat.var = 'Y.CENTER', size = 5, remove.na = FALSE) # show all grid names 
p14 <- p6 %>% add_label(dat = gridded.map.mean@data, lab.var = 'NAME', lon.var = 'X.CENTER', 
                        lat.var = 'Y.CENTER', size = 5, remove.na = TRUE) # only show the names of grid covered by field sampling 
cowplot::plot_grid(p13, p14, align = 'hv', ncol = 2, label_size = 28, labels = c("a", "b")) %>% suppressWarnings()
No description has been provided for this image
In [24]:
# Add the names of naive polygons and sample number in each naive polygon onto the map 
options(repr.plot.width = 15.4, repr.plot.height = 8.02 * 3)
p15 <- p4 %>% 
    add_label(dat = common.map.mean@data, lab.var = 'NAME', lon.var = 'X.CENTER', lat.var = 'Y.CENTER', size = 5) %>% 
    add_label(dat = common.map.mean@data, lab.var = 'sample.num', lon.var = 'X.CENTER', lat.var = 'Y.CENTER', size = 5, vjust = 3) %>% suppressWarnings()

# Add the sample number of each grid onto the map
p16 <- p7 %>% add_label(dat = gridded.map.mean@data, lab.var = 'sample.num', lon.var = 'X.CENTER', lat.var = 'Y.CENTER', size = 5, remove.na = F)

# Add the standard errors of shannon index in each grid on the map 
# Grids without labels mean that the samples are not enough for the calculation of standard error (at least two samples are required)
p17 <- p7 %>% add_label(dat = gridded.map.mean@data, lab.var = 'shannon_se', lon.var = 'X.CENTER', lat.var = 'Y.CENTER', size = 5, remove.na = F) 

# Display the figure
cowplot::plot_grid(p15, p16, p17, align = 'hv', ncol = 1, label_size = 28, labels = c("a", "b", "c")) %>% suppressWarnings()
No description has been provided for this image

Add points onto a map¶

In some cases, we want to display additional data onto the map processed by grid_map(). For example, we need to display both the mean and standard error of the Shannon index in each grid at the same time. The first option is to represent the mean by the color of the grid, and then use the add_label() function to directly add the value of the standard error to the grid. Another option is to use add_point() function, which adds points to each grid, with the size of the point representing the standard error.

In [25]:
# Add points onto a gridded map by using the standard error of shannon index in each grid
# Grids without point mean that the samples are not enough for the calculation of standard error (at least two samples are required)
options(repr.plot.width = 15.4, repr.plot.height = 8.02)
(p18 <- p7 %>% add_point(dat = gridded.map.mean@data, lab.var = 'shannon_se', 
                              lon.var = "X.CENTER", lat.var = "Y.CENTER", remove.na = TRUE) )# Remove all NA
No description has been provided for this image

You also can use the add_point() to visualize a numeric variable onto maps. For example, we can visualize the soil pH of each sampling site onto maps. The below is a simple example.

In [26]:
# Create a microgeo dataset. See subsequent tutorial for details. 
data(qtp)
map <- read_aliyun_map(adcode = c(540000, 630000, 510000))
dataset.dts <- create_dataset(mat = qtp$asv, ant = qtp$tax, met = qtp$met, map = map, 
                              phy = qtp$tre, env = qtp$env, lon = "longitude", lat = "latitude")
dataset.dts %<>% tidy_dataset()

ℹ [2024-01-04 17:10:15.457001] INFO ==> all samples fall within the map area!

ℹ [2024-01-04 17:10:15.464771] INFO ==> dataset has been created successfully!

ℹ [2024-01-04 17:10:15.470771] INFO ==> use `object %>% show_dataset()` to check the summary of dataset.
In [27]:
# Visualize the soil pH of each sampling site onto maps
options(repr.plot.width = 15.4 * 2, repr.plot.height = 8.02 * 3)
dat <- cbind(dataset.dts$met, dataset.dts$env)
p19 <- p1  %>% add_point(dat = dat, lab.var = 'pH') 
p20 <- p1  %>% add_point(dat = dat, lab.var = 'pH', fill = 'lightblue') 
p21 <- p1  %>% add_point(dat = dat, lab.var = 'pH', fill = rev(colorRampPalette(RColorBrewer::brewer.pal(11, 'RdYlGn'))(100)))
p22 <- p7  %>% add_point(dat = dat, lab.var = 'pH') 
p23 <- p9  %>% add_point(dat = dat, lab.var = 'pH')
p24 <- p11 %>% add_point(dat = dat, lab.var = 'pH')
cowplot::plot_grid(p19, p20, p21, p22, p23, p24, align = 'hv', ncol = 2, label_size = 28, 
                   labels = c("a", "b", "c", "d", "e", "f")) %>% suppressWarnings()
No description has been provided for this image

Add SpatRasters onto a map¶

In some cases, we need to visualize SpatRaster data. For example, spatial datasets (such as aridity index, annual average temperature, and soil pH) downloaded from public databases, and the results returned by some spatial interpolation functions and machine learning prediction functions in the microgeo R package. Regarding to the visualization of SpatRaster data, we have designed a function of add_spatraster(). Here are a few simple examples.

In [28]:
# Create a microgeo dataset. See subsequent tutorial for details. 
data(qtp)
map <- read_aliyun_map(adcode = c(540000, 630000, 510000))
dataset.dts <- create_dataset(mat = qtp$asv, ant = qtp$tax, met = qtp$met, map = map,
                              phy = qtp$tre, env = qtp$env, lon = "longitude", lat = "latitude")

ℹ [2024-01-04 17:10:23.413691] INFO ==> all samples fall within the map area!

ℹ [2024-01-04 17:10:23.421172] INFO ==> dataset has been created successfully!

ℹ [2024-01-04 17:10:23.427492] INFO ==> use `object %>% show_dataset()` to check the summary of dataset.
In [29]:
# Download phh2o from soilgrid database.
dataset.dts %<>% get_soilgrid(depth = 5, measures = c('phh2o'), out.dir = "../dev/dat/rundata")

✔ [2024-01-04 17:10:24.499287] SAVE ==> results have been saved to: object$spa$rast$his$phh2o
In [30]:
# Add the SpatRaster data of soil pH onto a map 
options(repr.plot.width = 15.4, repr.plot.height = 8.02)
(p25 <- p1 %>% add_spatraster(spat.raster = dataset.dts$spa$rast$his$phh2o))
No description has been provided for this image

Notably, the add_spatraster() function would fill the whole study area. However, in some cases, we only focus one or two ecosystems, for example, grassland, forest or barren. To address this problem, we also have implemented a function of mask_spatraster_by_cla(), which can mask off those areas that we are not interested in. Here are a few simple examples.

In [31]:
# Download classification MODIS metrics of research region. 
# Please replace the <username> and <passwd> with real value
dataset.dts %<>% get_modis_cla_metrics(username = "username", password = "passwd", out.dir = "../dev/dat/rundata")

ℹ [2024-01-04 17:10:27.520267] INFO ==> preparing MODIS product list for searching...

ℹ [2024-01-04 17:10:27.535475] INFO ==> searching avaliable MODIS products...

ℹ [2024-01-04 17:10:27.542682] INFO ==> current product (1/1): MCD12Q1 (LC_Type1--> 2022-01-01 to 2022-12-31)

ℹ [2024-01-04 17:10:29.238166] INFO ==> find 8 files with 0.09 GB in total...

ℹ [2024-01-04 17:10:29.247663] INFO ==> downloading all avaliable MODIS products[skip if the file exists]...

ℹ [2024-01-04 17:10:29.257661] INFO ==> preparing the PTVs (Product, Time, Version) for merging remote-sensing images...

ℹ [2024-01-04 17:10:29.28531] INFO ==> converting hdf files to tif files...

ℹ [2024-01-04 17:10:29.293728] INFO ==> current product (1/1): MCD12Q1 (convert 8 hdf files into 1 tif files using 1 threads)

ℹ [2024-01-04 17:10:29.506624] INFO ==> collecting all merged image files...

ℹ [2024-01-04 17:10:29.515717] INFO ==> current measure (1/1): LC_Type1_061

ℹ [2024-01-04 17:10:29.588037] INFO ==> reprojecting the CRS of SpatRaster to epsg:+proj=longlat +datum=WGS84 +no_defs, it takes a while...

✔ [2024-01-04 17:11:01.754878] SAVE ==> results have been saved to: object$spa$rast$cla$LC_Type1
In [32]:
# Mask the results by using grassland(10); See `?mask_spatraster_by_cla()` for details about the classification codes of ecosystems
options(repr.plot.width = 15.4, repr.plot.height = 8.02)
m1 <- mask_spatraster_by_cla(tar.spat = dataset.dts$spa$rast$his$phh2o, ref.spat = dataset.dts$spa$rast$cla$LC_Type1, use.class = 10)
(p26 <- p1 %>% add_spatraster(spat.raster = m1))
No description has been provided for this image
In [33]:
# Mask the results by using barren(16); See `?mask_spatraster_by_cla()` for details about the classification codes of ecosystems
options(repr.plot.width = 15.4, repr.plot.height = 8.02)
m2 <- mask_spatraster_by_cla(tar.spat = dataset.dts$spa$rast$his$phh2o, ref.spat = dataset.dts$spa$rast$cla$LC_Type1, use.class = 16)
(p27 <- p1 %>% add_spatraster(spat.raster = m2))
No description has been provided for this image
In [34]:
# Mask the results by using forest (1,2,3,4,5); See `?mask_spatraster_by_cla()` for details about the classification codes of ecosystems
options(repr.plot.width = 15.4, repr.plot.height = 8.02)
m3 <- mask_spatraster_by_cla(tar.spat = dataset.dts$spa$rast$his$phh2o, ref.spat = dataset.dts$spa$rast$cla$LC_Type1, use.class = c(1,2,3,4,5))
(p28 <- p1 %>% add_spatraster(spat.raster = m3))
No description has been provided for this image
In [35]:
# Mask the results by using grassland (10) and barren (16); See `?mask_spatraster_by_cla()` for details about the classification codes of ecosystems
options(repr.plot.width = 15.4, repr.plot.height = 8.02)
m4 <- mask_spatraster_by_cla(tar.spat = dataset.dts$spa$rast$his$phh2o, ref.spat = dataset.dts$spa$rast$cla$LC_Type1, use.class = c(10, 16))
(p29 <- p1 %>% add_spatraster(spat.raster = m4))
No description has been provided for this image

Add contours onto a map¶

You also can convert the SpatRaster returned by spatial interpolation functions or machine learning prediction functions to contours, and further add these contours onto maps by using add_contour(). Here are a few simple examples.

In [36]:
# Create a microgeo dataset. See subsequent tutorial for details. 
data(qtp)
map <- read_aliyun_map(adcode = c(540000, 630000, 510000))
dataset.dts <- create_dataset(mat = qtp$asv, ant = qtp$tax, met = qtp$met, map = map,
                              phy = qtp$tre, env = qtp$env, lon = "longitude", lat = "latitude") 
dataset.dts %<>% rarefy_count_table() 
dataset.dts %<>% tidy_dataset()
dataset.dts %<>% calc_alpha_div(measures = c("observed", "shannon"))

ℹ [2024-01-04 17:11:48.354824] INFO ==> all samples fall within the map area!

ℹ [2024-01-04 17:11:48.364121] INFO ==> dataset has been created successfully!

ℹ [2024-01-04 17:11:48.372576] INFO ==> use `object %>% show_dataset()` to check the summary of dataset.

ℹ [2024-01-04 17:11:50.699992] INFO ==> the ASV/gene abundance table has been rarefied with a sub-sample depth of 5310

✔ [2024-01-04 17:11:53.637593] SAVE ==> results have been saved to: object$div$alpha
In [37]:
# Perform kriging interpolation for shannon diversity indices.
kri.rst.shannon <- interp_kri(map = dataset.dts$map, met = dataset.dts$met, dat = dataset.dts$div$alpha, 
                              var = 'shannon', model = 'Mat', trim.dup = TRUE)

! [2024-01-04 17:11:53.963551] WARN ==> only use 420 out of 1244 sampling sites for interpolation!


[using ordinary kriging]
In [38]:
# Add contours on the map 
options(repr.plot.width = 15.4 * 2, repr.plot.height = 8.02)
p30 <- p1 %>% add_contour(spat.raster = kri.rst.shannon, nlevels = 100)
p31 <- p1 %>% add_spatraster(spat.raster = kri.rst.shannon) %>% add_contour(spat.raster = kri.rst.shannon, nlevels = 10, color = 'white', show.labels = TRUE)
cowplot::plot_grid(p30, p31, align = 'hv', ncol = 2, label_size = 28, labels = c("a", "b")) %>% suppressWarnings()
No description has been provided for this image

Add north arrow onto a map¶

We have designed a function of add_north_arrow() to add a north arrow onto the map. You can use ?add_north_arrow() to view the built-in document. Here are a few simple examples.

In [39]:
# Add a north arrow on the map
options(repr.plot.width = 15.4 * 2, repr.plot.height = 8.02)
p32 <- p29 %>% add_north_arrow() %>% suppressWarnings()
p33 <- p30 %>% add_north_arrow() %>% suppressWarnings()
cowplot::plot_grid(p32, p33, align = 'hv', ncol = 2, labels = c("a", "b"), label_size = 28) %>% suppressWarnings()
No description has been provided for this image

Add spatial-aware scale bar onto a map¶

We have designed a function of add_scale_bar() to add a spatial-aware scale bar onto the map. You can use ?add_scale_bar() to view the built-in document. Here are a few simple examples.

In [40]:
# Add a spatial-aware scale bar on the map
options(repr.plot.width = 15.4 * 2, repr.plot.height = 8.02)
p34 <- p32 %>% add_scale_bar()
p35 <- p33 %>% add_scale_bar()
cowplot::plot_grid(p34, p35, align = 'hv', ncol = 2, labels = c("a", "b"), label_size = 28) %>% suppressWarnings() %>% suppressMessages()
No description has been provided for this image

Add CRS onto a map before final visualization¶

We have designed a function of add_crs() to add a coordinate reference system to the map before final visualization. You can use ?add_crs() to view the built-in document. Here are a few simple examples.

In [41]:
# Add a coordinate reference system to the map before final visualization
options(repr.plot.width = 13.4 * 2, repr.plot.height = 8.02)
p36 <- p34 %>% add_crs()
p37 <- p35 %>% add_crs()
cowplot::plot_grid(p36, p37, align = 'hv', ncol = 2, labels = c("a", "b"), label_size = 28) %>% suppressWarnings() %>% suppressMessages()
No description has been provided for this image

Add sampling sites onto a map¶

In the large-scale microbial biogeographical studies, another application of microgeo R package is to draw sample point distribution maps. In the current version, we provide a flexible function of add_sampl_site(), which can render sampling points on the map. Here are a few simple examples.

In [42]:
# Create a microgeo dataset. See subsequent tutorial for more details.
data(qtp)
map <- read_aliyun_map(adcode = c(540000, 630000, 510000))
dataset.dts <- create_dataset(mat = qtp$asv, ant = qtp$tax, met = qtp$met, map = map,
                              phy = qtp$tre, env = qtp$env, lon = "longitude", lat = "latitude") 
dataset.dts %<>% get_ai(out.dir = "../dev/dat/rundata") 
dataset.dts %<>% get_his_bioc(res = 10, out.dir = "../dev/dat/rundata") 
dataset.dts %<>% get_soilgrid(depth = 5, measures = c('phh2o'), out.dir = "../dev/dat/rundata") 
dataset.dts %<>% extract_data_from_spatraster(type = 'his') 
dataset.dts %<>% tidy_dataset()

ℹ [2024-01-04 17:12:08.654725] INFO ==> all samples fall within the map area!

ℹ [2024-01-04 17:12:08.664679] INFO ==> dataset has been created successfully!

ℹ [2024-01-04 17:12:08.674071] INFO ==> use `object %>% show_dataset()` to check the summary of dataset.

✔ [2024-01-04 17:12:12.710911] SAVE ==> results have been saved to: object$spa$rast$his$AI

✔ [2024-01-04 17:12:20.137347] SAVE ==> results have been saved to: object$spa$rast$his(19 variables)

✔ [2024-01-04 17:12:26.392774] SAVE ==> results have been saved to: object$spa$rast$his$phh2o

✔ [2024-01-04 17:12:26.411119] SAVE ==> results have been saved to: object$spa$tabs
In [43]:
# Add two grouping variable in `dataset.dts$met`
dataset.dts$met$GA <- cut(dataset.dts$spa$tabs$AI, breaks = c(-Inf, 0.2, 0.5, Inf), labels = c("Arid", "Semi-arid", "Humid"))
dataset.dts$met$GM <- cut(dataset.dts$spa$tabs$Bio1, breaks = c(-Inf, -1, 1, Inf), labels = c("Cold", "Moderate", "Warm"))
In [44]:
# Add sampling sites without colors 
options(repr.plot.width = 13.4 * 2, repr.plot.height = 8.02 * 3)
s1 <- p4 %>% add_sampl_site(met = dataset.dts$met, color.val = "gray30", fill.val = "white", point.alpha = 0.5) %>% 
    add_scale_bar() %>% add_north_arrow() %>% add_crs()
s2 <- p7 %>% add_sampl_site(met = dataset.dts$met, color.val = "gray30", fill.val = "white", point.alpha = 0.5) %>% 
    add_scale_bar() %>% add_north_arrow() %>% add_crs()
s3 <- p11 %>% add_sampl_site(met = dataset.dts$met, color.val = "gray30", fill.val = "white", point.alpha = 0.5) %>% 
    add_scale_bar() %>% add_north_arrow() %>% add_crs()
s4 <- p9 %>% add_sampl_site(met = dataset.dts$met, color.val = "gray30", fill.val = "white", point.alpha = 0.5) %>% 
    add_scale_bar() %>% add_north_arrow() %>% add_crs()
s5 <- p25 %>% add_sampl_site(met = dataset.dts$met, color.val = "gray30", fill.val = "white", point.alpha = 0.5) %>% 
    add_scale_bar() %>% add_north_arrow() %>% add_crs()
s6 <- p29 %>% add_sampl_site(met = dataset.dts$met, color.val = "gray30", fill.val = "white", point.alpha = 0.5) %>% 
    suppressWarnings() %>% add_scale_bar() %>% add_north_arrow() %>% add_crs()
cowplot::plot_grid(s1, s2, s3, s4, s5, s6, align = 'hv', ncol = 2, labels = c("a", "b", "c", "d", "e", "f"), label_size = 28) %>% suppressWarnings()
No description has been provided for this image
In [45]:
# Add sampling sites with color and shapes
options(repr.plot.width = 13.4 * 2, repr.plot.height = 8.02)
s7 <- p1 %>% add_sampl_site(met = dataset.dts$met, color.var = "GA", color.val = "gray30",point.alpha = 0.5) %>%  
    add_scale_bar() %>% add_north_arrow() %>% add_crs()
s8 <- p1 %>% add_sampl_site(met = dataset.dts$met, color.var = "GA", shape.var = "GA", color.val = "gray30", point.alpha = 0.5) %>%
    add_scale_bar() %>% add_north_arrow() %>% add_crs()
cowplot::plot_grid(s7, s8, align = 'hv', ncol = 2, labels = c("a", "b"), label_size = 28) %>% suppressWarnings()
No description has been provided for this image
In [46]:
# Facet by one grouping variable
options(repr.plot.width = 13.43 * 2, repr.plot.height = 7.9 * 3)
s9 <- p1 %>% add_sampl_site(met = dataset.dts$met, color.var = "GA", shape.var = "GA", color.val = "gray30", point.alpha = 0.5, facet.by.color = T, facet.col.nums = 1) %>%
      add_scale_bar() %>% add_north_arrow() %>% add_crs()
s10 <- p29 %>% add_sampl_site(met = dataset.dts$met, color.var = "GA", shape.var = "GA", color.val = "gray30", point.alpha = 0.5, facet.by.color = T, facet.col.nums = 1) %>%
       add_scale_bar() %>% add_north_arrow() %>% add_crs()
cowplot::plot_grid(s9, s10, align = 'hv', ncol = 2, labels = c("a", "b"), label_size = 28) %>% suppressWarnings()
No description has been provided for this image
In [47]:
# Facet by two grouping variables [on basic map]
options(repr.plot.width = 13.43 * 3, repr.plot.height = 7.23 * 3)
p1 %>% add_sampl_site(met = dataset.dts$met, color.var = "GA", shape.var = "GM", color.val = "gray30", point.alpha = 0.5, facet.by.color = T, facet.by.shape = T) %>%
       add_scale_bar() %>% add_north_arrow() %>% add_crs()
No description has been provided for this image
In [48]:
# Facet by two grouping variables [on SpatRaster]
options(repr.plot.width = 13.43 * 3, repr.plot.height = 7.9 * 3)
p29 %>% add_sampl_site(met = dataset.dts$met, color.var = "GA", shape.var = "GM", color.val = "gray30", point.alpha = 0.5, facet.by.color = T, facet.by.shape = T) %>%
        add_scale_bar() %>% add_north_arrow() %>% add_crs()
No description has been provided for this image